######## INFO ########

# PROJECT
# Replication files for: Economic Impact of the South China Sea Dispute in China-Philippines
#                        Relations from 2012 to 2016: A DM-LFM analysis  

# Authors: Dianyi Yang
# R Script
# Purpose: This script produces the remaining figures for exports to China.
# Created: 22 Oct 2023
# Updated: 29 Nov 2023
# Inputs: stored/dmlfm.rdata
#         custom_functions/DM-LFM.r
#         stored/dmlfm_placebo.rdata
#         stored/leave-one-out.rdata
# Outputs: output/dmlfm_trend.pdf
#          output/dmlfm_trend.png
#          output/dmlfm_effect.pdf
#          output/dmlfm_effect.png
#          output/placebo_trend.pdf
#          output/placebo_trend.png
#          output/loo_trend.pdf
#          output/loo_trend.png
#          output/loo_effect.pdf
#          output/loo_effect.png


######## SETUP ########
rm(list = ls()) # clear workspace

need <- c("Synth", "reshape2","devtools", "synthdid",'tidyverse',
          'bpCausal', "fixest","modelsummary",
          'kableExtra') # list packages needed
have <- need %in% rownames(installed.packages()) # checks packages you have
if(any(!have)) install.packages(need[!have]) # install missing packages
#devtools::install_github("synth-inference/synthdid") #manually install sdid if the codes above fail to do so automatically
invisible(lapply(need, library, character.only=T)) # load needed packages

setwd(gsub('/code','',dirname(rstudioapi::getSourceEditorContext()$path))) #set wd to root directory
list.files()
##################DM-LFM####################################
load('stored/dmlfm.rdata')
source('custom_functions/DM-LFM.r')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
#plot.DMLFM(out, ci.level=0.90)
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI
class(result.use90) <- 'DMLFM'

#trend plot
ggplot(result.use$est.eff, aes(x=time))+
  geom_ribbon(aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "95% CI"), alpha = 0.5)+
  geom_ribbon(data=result.use90$est.eff,aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "90% CI"), alpha = 0.5)+
  geom_line(aes(y=observed, color='Observed'),linewidth=0.75,alpha=0.5)+
  geom_line(aes(y=estimated_counterfactual, color='Counterfactual'),linewidth=0.75,alpha=0.5)+
  geom_vline(xintercept=0)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+ 
  scale_color_manual(name= '',
                     breaks=c('Observed','Counterfactual'),
                     values=c('Observed'="#00BFC4", "Counterfactual"="#F8766D"))+
  labs(x="Months since treatment", y='Log(Filipino exports to China)')+
  scale_fill_manual(name=c("", ''),
                    breaks=c('95% CI', '90% CI'),
                    values=c('lightgrey', 'grey'))  
ggsave('output/dmlfm_trend.pdf', width = 10, height=5) #Figure 3
ggsave('output/dmlfm_trend.png', width = 10, height=5)

#effect plot
ggplot(result.use$est.eff, aes(x=time, y=estimated_ATT))+
  geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u, fill = "95% CI"), alpha = 0.5)+
  geom_ribbon(data=result.use90$est.eff,aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u, fill = "90% CI"), alpha = 0.5)+
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +
  theme_bw()+labs(x="Months since treatment",y="ATT with 90% and 95% CIs") + geom_line(linewidth=0.75) +
  scale_fill_manual(name=c("", ''),
                    breaks=c('95% CI', '90% CI'),
                    values=c('lightgrey', 'grey'))         #effect plot
ggsave('output/dmlfm_effect.pdf', width = 10, height=5)  #Figure 5
ggsave('output/dmlfm_effect.png', width = 10, height=5)

#Placebo test (DMLFM)
load('stored/dmlfm_placebo.rdata')
result.use <- effSummary(out) #extract effect information
class(result.use) <- 'DMLFM'
result.use90 <- effSummary(out, ci.level = 0.90) #for 90% CI

#trend plot
ggplot(result.use$est.eff, aes(x=time-25))+
  geom_ribbon(aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "95% CI"), alpha = 0.5)+
  geom_ribbon(data=result.use90$est.eff,aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u, fill = "90% CI"), alpha = 0.5)+
  geom_line(aes(y=observed, color='Observed'),linewidth=0.75,alpha=0.5)+
  geom_line(aes(y=estimated_counterfactual, color='Counterfactual'),linewidth=0.75,alpha=0.5)+
  geom_vline(xintercept=-25)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+ 
  scale_color_manual(name= '',
                     breaks=c('Observed','Counterfactual'),
                     values=c('Observed'="#00BFC4", "Counterfactual"="#F8766D"))+
  labs(x="Months since actual treatment", y='Log(Filipino exports to China)')+
  scale_fill_manual(name=c("", ''),
                    breaks=c('95% CI', '90% CI'),
                    values=c('lightgrey', 'grey'))  
ggsave('output/placebo_trend.pdf', width = 10, height=5) #Figure 7
ggsave('output/placebo_trend.png', width = 10, height=5)

#effect plot (90% and 95% CIs)
ggplot(result.use$est.eff, aes(x=time-25, y=estimated_ATT))+
  geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), fill = "lightgrey", alpha = 0.5)+
  geom_ribbon(data=result.use90$est.eff,aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), fill = "grey", alpha = 0.5)+
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey") +
  geom_vline(xintercept = -25, linetype = "dashed", color = "grey") +
  geom_line(linewidth=0.75)+
  theme_bw()+labs(x="Months since treatment",y="ATT with 90% and 95% CIs")

#leave-one-out test
load('stored/leave-one-out.rdata')

narrow <- plot %>%
  group_by(time) %>%
  summarize(counterfactual_ci_l = max(counterfactual_ci_l), counterfactual_ci_u = min(counterfactual_ci_u))

#trend
ggplot(plot, aes(x=time))+
  geom_ribbon(aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u), alpha = 0.5, fill='lightgrey')+
  geom_ribbon(data=narrow,aes(ymin = counterfactual_ci_l, ymax = counterfactual_ci_u), alpha = 0.5, fill='grey')+
  geom_line(data=subset(plot, miss_unit!='All'),aes(y=estimated_counterfactual, color= miss_unit),linewidth=0.5,alpha=0.5)+
  geom_line(data=subset(plot, miss_unit=='All'),aes(y=estimated_counterfactual, color='Counterfactual (full pool)'),linewidth=1,alpha=1)+
  geom_line(aes(y=observed, color='Observed'),linewidth=1,alpha=1)+
  geom_vline(xintercept=0)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+
  labs(x="Months since treatment", y='Log(Filipino exports to China)')+
  scale_color_discrete(name='Key',
                       breaks=c('Observed', 'Counterfactual (full pool)', unique(plot$miss_unit)))
ggsave('output/loo_trend.pdf', width = 10, height=5) #Figure 9
ggsave('output/loo_trend.png', width = 10, height=5)

#effect
narrow <- plot %>%
  group_by(time) %>%
  summarize(estimated_ATT_ci_l = max(estimated_ATT_ci_l), estimated_ATT_ci_u = min(estimated_ATT_ci_u))

ggplot(plot, aes(x=time))+
  geom_ribbon(aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), alpha = 0.5, fill='lightgrey')+
  geom_ribbon(data=narrow,aes(ymin = estimated_ATT_ci_l, ymax = estimated_ATT_ci_u), alpha = 0.5, fill='grey')+
  geom_line(data=subset(plot, miss_unit!='All'),aes(y=estimated_ATT, color= miss_unit),linewidth=0.5,alpha=0.5)+
  geom_line(data=subset(plot, miss_unit=='All'),aes(y=estimated_ATT, color='ATT with full pool'),linewidth=1,alpha=1)+
  geom_vline(xintercept=0)+theme_minimal() + theme(panel.background = element_rect(fill = "white"))+
  labs(x="Months since treatment", y='ATT with widest and narrowest 95% CIs')+
  scale_color_discrete(name='Key',
                       breaks=c('ATT with full pool', unique(plot$miss_unit)))
ggsave('output/loo_effect.pdf', width = 10, height=5) #Figure 11
ggsave('output/loo_effect.png', width = 10, height=5)